The competition between the formation and destruction of 
coronene clusters under interstellar conditions is investigated 
theoretically. The unimolecular nucleation of neutral clusters is 
simulated with an atomic model combining an explicit classical 
force field and a quantum tight-binding approach. Evaporation 
rates are calculated in the framework of the phase space the- 
ory and are inserted in an infrared emission model and com- 
pared with the growth rate constants. It is found that, in inter- 
stellar conditions, most collisions lead to cluster growth. The 
time evolution of small clusters (containing up to 312 carbon 
atoms) was specifically investigated under the physical condi- 
tions of the northern photodissociation region of NGC 7023. 
These clusters are found to be thermally photoevaporated much 

- , faster than they are reformed, thus providing an interpretation 
for the lowest limit of the interstellar cluster size distribution 
, inferred from observations. The effects of ionizing the clusters 

PsJ ' and density heterogeneities are also considered. Based on our 
results, the possibility that PAH clusters could be formed in 
(]J PDRs is critically discussed. 
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Abstract. 



1. Introduction 

Many astronomical objects show a distinct set of emission 
bands in the mid-infrared range, known collectively as the 
unidentified infrared emission features. The most intense of 
these features, the Aromatic Infrared Bands (AIBs) falling 
at 3.3, 6.2, "7.7", 8.6, 11.3 and 12.7 jum, are the signa- 
tures of aromatic CC and CH bonds. They are observed sys- 
tematically from different regions of the interstellar medium 
(ISM) irradiated by UV photons. The carriers of these bands 
have been identified as Polycych c Aromatic Hydrocar bons 
(PAHs) some twenty years ago by iLeger & Pugetl ( 1 19841) and 

[Allamandola et al.. (.1985.). 

iBoulanger et alJ (ll99Cl) and lBernard et al J lll993h analysed 
the emission in the IRAS photometric bands measured for sev- 
eral molecular clouds, and suggested that free-flying PAHs 
are produced by photoevaporation of larger grains. The ob- 
servation of a mid-infrared cont inuum in the re flection nebula 
Ced 201 has been attributed bv iCesarskv et alJ (.2000.) to Very 
Small carbonaceous Grains (VSGs) producing the AIB carri- 
ers. More recently, in their study of the two reflection nebu- 
lae NGC 7023 and p Ophiucus-SR3. lRapacioh et alJ (l2005bl) 
showed that PAH molecules are produced by the destruction 
of small carbonaceous grains inside molecular clouds, these 
grains being interpreted as PAH clusters. A minimal size of 
400 carbon atoms per cluster was inferred from the analysis 
of astronomical observations (Rapacioli et al. 2005b). A plau- 
sible scenario is that these grains undergo photo-evaporation at 
the edge of clouds, thus leading to isolated PAHs. These ob- 
servations motivate theoretical and experimental studies on the 
physical properties of PAH clusters, such as structure, stability 
or aggregation. 
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Small clusters of PAH molecules have been studied by sev- 
eral a uthors in the past , both experimentally (Benharash etal J 



1999'. ISong et alJl2003l: IPiuzzi et alJl2002l iMiUer et alJ 


19841 


and theoretically ("van de Waal' 


1983"; 'Gonzales & Lim| 


199?; 


Marzec 2000; 


Piuzzi et al. 2002; Grimme 2004). Most of these 



studies have focused on rather small species, usually con- 
taining no more than a few tens of carbon atoms, far be- 
low the limit of 400 ca rbon atoms per VSG as suggested by 
iRapacioli et alJ ll2005bl) . Studies on larger systems have been 
initiated only lately, motivated by the astrophysical context. 
Clusters of coronene molecules containing u p to 20 units have 
been produced in a gas ag gregation source JBrechignac et alJ 
I2005t1schmidt et alJl2006 ). Interestingly, the authors showed 
that the dissociation of coronene clusters induced by multi- 
photon excitation of a near-UV laser proceeds via the ejection 
of van der Waals-bonded coronene units. 

In a previous work jRapacioli et alJ Eo05al) . we investi- 
gated the stable structures of PAH clusters containing up to 800 
carbon atoms (32 coronene molecules). The considered PAHs 
ranged from pyrene (CieHio) to circumcoronene (C54H18). The 
most stable cluster structures generally result from stacking the 
PAH molecules, first in a one-dimensional pattern at small clus- 
ter sizes. Above some critical size (depending on the PAH it- 
self), structures more stable than the single stack are found, 
consisting of shorter stacks lying next to each other, and grow- 
ing as two-, eventually three-dimensional structures. In large 
PAH clusters, the multiple ways of arranging the short stacks 
lead to a competition between stable structures. 

Our present interest lies in the formation and destruction 
of PAH clusters in the ISM, and particularly the competition 
between these two processes. The coronene molecule (C24H12) 
provides a convenient example of large PAHs, especially from 
the computational point of view, and is used in the present work 
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as a prototype for interstellar PAHs. We wish to address here 
the following questions; 

(i) What are the astrophysical conditions that favour the nucle- 
ation of PAH clusters? 

(ii) How much energy is required for thermal evaporation of 
these clusters? 

(iii) How do the rates of formation and destruction compare for 
these clusters under interstellar conditions? 

In this paper, we focus on small neutral coronene clusters 
containing up to 312 carbon atoms, or 13 molecules. In particu- 
lar, we attempt to offer a physical explanation for the existence 
of the observed lower limit to the interstellar PAH cluster size 
distribution of 400 carbon atoms per cluster iRaoacioli et alJ 
l2005bh . 

Addressing these questions requires a variety of theoret- 
ical methods to be used. Nucleation processes are explicitly 
simulated with an extended version of the tight-bin ding (TB) 
model initially developed bv lVan Oanh et al.1 (l2002h for single 
PAHs. Concerning the thermal stability of coronene clusters, 
the rates of unimolecular evaporation as well as the photoemis- 
sion rates in the infrared are obtained and compared. The evap- 
oration rates are obtained in the statistical framework of the 
phase space theory (Chesnavich & Bowers 1977). The compe- 
tition between evaporation and IR emission is calculated using 
a master equation modeling based on microcanonical statistics 
( Joblin et al. 2002). In the light of these calculations, we dis- 
cuss several chemical scenarios to interpret th e observations of 
the reflection nebulae jRapacioli et alJ2005bl) . 

The paper is organized as follows. In the next section, we 
describe the atomic model and discuss the numerical simula- 
tion of coronene clusters growth by unimolecular nucleation. 
In section|3 the evaporation rates of coronene clusters are cal- 
culated. In section 0] we study the competition between evap- 
oration and IR emission. In section |5] the photophysics of a 
coronene cluster in a radiation field is simulated. In section|6l 
we study the competition between formation and destruction of 
coronene clusters under interstellar conditions. We also high- 
light the astrophysical implications of our results for possible 
models of chemical evolution in the photodissociation regions. 
In particular, the possible roles of ionization and density het- 
erogeneities in the PDR in the existence and stability of PAH 
clusters are discussed. Some conluding remarks are given in 
Sec.E] 

2. Cluster formation process 

In this section, we consider the formation of neutral coronene 
clusters by collision of a single molecule with a smaller cluster 
or with another single molecule. In the ISM, collisions occur in 
a very rarefied gas and the time lapse between successive UV 
absorption events is in general much longer than the IR emis- 
sion timescale. Therefore the PAH molecules and their clusters 
have enough time to cool down through IR emission, and can 
be considered as initially cold in their ground electronic and vi- 
brational states. This energy will be transferred to intermolec- 
ular modes, then to the intramolecular modes. Intramolecular 
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vibrational energy redistribution (IVR) is likely to be the ma- 
jor source for dissipating the initial energy, having thus strong 
consequences on the collision products. Hence it is important to 
take this redistribution into account in molecular simulations. 

In this paper, we will consider that a new cluster is formed 
upon collision if it has time to reach thermal equilibrium be- 
fore further dissociating. The relaxation of a thermalised clus- 
ter heated by collision or by photoabsorption is discussed in 
section|3] We describe hereafter the atomic model that was im- 
plemented for the present problem. 

2.1. Atomic model 

Describing the transfer of the collisional energy toward inter- 
molecular and intramolecular modes requires both types of de- 
grees of freedom to be correctly accounted for Clusters of 
coronene molecules are modeled using the combination of an 
explicit classical force field to describe the intermolecular in- 
teractions, along with a quantum tight-binding model for the 
intramolecular interactions within each molecule. The total 
Hamiltonian for a system with A^^oi PAH molecules is written: 

^=I]^La+^^inte.- (1) 

In this equation H'. , is the intramolecular Hamiltonian of 

^ intra 

molecule /, which only explicitly depends on the positions of 
the atoms in this m olecule. We use for H[ the tight-binding 
model developed bv lVan Oanh et alJ (120021) . This simple quan- 
tum model was parameterized to reproduce energetic and vi- 
brational properties of isolated PAH molecules. We did not 
change its original parameters. 

//inter IS the intcrmolccular potential resulting from the 
repulsion-dispersion forces as well as the electrostatic in- 
teractions between partial charges on the a toms of differ- 
ent molecules. Following our previous work ("Rapa cioli et alJ 
,2005a) . we use a simple Lennard-Jones (LJ) -i- point charges 
potential, the charges being determined to reproduce Density 
Functional Theory (DFT) calculations on an isolated molecule. 
The LJ parameters are taken from van de Waal (1983). Further 
details about the intermolecular potential, including a discus- 
sion about the sensitivity of the results with respect to the po- 
tential parameters, can be found in .Rapacioli et al.. (,2005a.) . 

2.2. Simulation details 

The classical equations of motion on the Born-Oppenheimer 
potential energy surface have been numerically integrated us- 
ing the fifth-order Adams-Moulton predictor-corrector algo- 
rithm. The time step was chosen between 0.05 fs and 0.15 fs 
depending on the collision energy. The computational effort re- 
quired to diagonalize the multiple TB Hamiltonians was dis- 
tributed among many processors on parallel computers in the 
OpenMP framework. As stated previously, the reactants are as- 
sumed to be cold, their vibrational energies being taken as less 
than 80 K (much lower than the collisional energy). The phys- 
ical parameters, which characterize a collision in the centre of 
mass reference frame, are the colhsional energy (translational 
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plus rotational components), the impact parameter, the angular 
momentum of each fragment and their initial orientations. The 
statistical information obtained from the simulations was lim- 
ited due to the still relatively heavy numerical cost involved in 
computing the TB energies and forces. Therefore we restricted 
our study to the following representative situations, assuming 
this restriction has a minor impact on the collision products: 

- all collisions are at zero impact parameter (no orbital angu- 
lar momentum); 

- the total angular momentum of the system (incident 
molecule-Hcluster) is zero; 

- the collisional energy is equally distributed between the 
translational energy and the rotational energy of the reac- 
tants. 

At a given collisional energy, the free mechanical parameters 
are thus reduced to the initial orientations of the molecule and 
the cluster and the direction of the angular momentum of the 
cluster We did not attempt to determine quantitative cross sec- 
tions for the collision processes. 

About 30 trajectories have been performed for each colU- 
sional energy in order to sample these parameters and to gather 
statistics on the outcome of the collisions. For simplicity the 
colUsion of reactants with sizes p and p' , respectively, leading 
to a set of fragments with sizes d, d' ,d" ... will be referred to as 
p + p' ^ d + d' + d" + ■ ■ ■, where the parents and products are 
ordered by increasing size. 

In practice the simulation ends when either one of the two 
following situations happens: 

- after the collision, fragmentation occurs and one or more 
molecules are ejected far from the others; 

- the instantaneous intra- and intermolecular temperatures 
cross each other, meaning that the product cluster has ab- 
sorbed the collisional energy and has reached thermal equi- 
librium. 

Other situations have not been observed. 
2.3. Results 

Figure ^ shows the nucleation probability of simulated colli- 
sions between two coronene molecules as a function of colli- 
sional energy, that is the chance that a collision leads to the 
formation of a dimer as 1 -H 1 -^2. The probability that the col- 
lision fails in producing a dimer (l-i-l — > l-i-l), also plotted 
in Fig.^ complements the nucleation probability. 

These results show that the probability of forming a dimer 
decreases with increasing collisional energy. This is not sur- 
prising, since collisions at higher energy tend to make the two 
molecules bounce of each other Nucleation is efficient, espe- 
cially at low energies. For instance, the probability of form- 
ing a coronene dimer exceeds 50% when the collisional energy 
is lower than 3.6 eV. This rather high energy threshold corre- 
sponds to the average collisional energy of a gas at a tempera- 
ture of 13500 K. 

The probabilities of successfully forming a larger cluster 
or of destroying the initial cluster from a collision between a 
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Fig. 1. Probabilities pi and /^JJ that a collision between two 
coronene molecules forms a dimer (1 -H 1 ^ 2) or not (1 -t- 1 — > 
1 -H 1), as a function of the collisional energy. Each data point 
is averaged over 30 simulations. Lines are drawn to guide the 
eye. 

single coronene and a stack of 2 or 3 molecules are shown in 
Figs.|5]and|3 respectively. For these clusters, the probability of 
a neutral event, that is when the collision products are identical 
to the parents in terms of their size, is also plotted. The mul- 
tifragmentation of the collision product into single molecules 
only was not observed in the case of the tetramer. 

As in the case of dimer formation, the probability of grow- 
ing the initial cluster decreases monotonically with increasing 
collisional energy. A monotonic increase is also observed for 
the probability of ending with smaller clusters. This allows us 
to define a critical energy for which these two probabilities are 
the same. This energy approximately equals 12.4 eV for the 
1-1-2 reaction and 14.2 eV for 1-1-3. The increase of the crit- 
ical energy with cluster size is expected since a larger number 
of modes are available to convert the collisional energy. In our 
simulations, the range of collisional energies was not chosen 
on the basis of astrophysical values but only to allow enough 
statistics to be gathered in the calculations. Even though the 
mean collisional energies encoutered in the interstellar medium 
are usually much lower (typically in the 0.01-0. 1 eV range), we 
extrapolate below the values corresponding to the smallest cal- 
culated collisional energies to collisions with lower energies. 
This approximation leads to some underestimation of the clus- 
ter growth, which we estimate to be no more than 25% since the 
sticking probability is at least 0.75 (as seen from Figs.^0. 

3. Evaporation of a neutral isolated cluster 

A thermalised isolated coronene cluster with an excitation en- 
ergy iitot can release a part of this excess energy via dissoci- 
ation, evaporation and IR emission. We will assume that the 
cluster remains on the ground state electronic potential surface. 
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Fig. 2. Probabilities pi, pi, and p_ that a collision of a single 
coronene on a dimer forms a trimer (1 + 2 — > 3), leaves a dimer 
(1+2 — > 1 + 2), or fragments into 3 molecules (1 + 2 — > 
1 + 1 + 1), as a function of the collisional energy. Each data 
point is averaged over 30 simulations. Lines are drawn to guide 
the eye. 



Fig. 3. Probabilities pi, pi, and pL that a collision of a single 
coronene on a trimer forms a tetramer (1 + 3 — > 4), leaves a 
trimer (1 + 3 — > 1+3), destroys the inital cluster (1+3 — > 2 + 2 
or 1+3 — > 1 + 1 +2), as a function of the collisional energy. Each 
data point is averaged over 30 simulations. Lines are drawn to 
guide the eye. 



In this section, we study the evaporation mechanism, which in 
turn requires an estimation of the characteristic time for the en- 
ergy redistribution between inter- and intramolecular modes. In 
the next section, we study the competition between the different 
relaxation mechanisms involving molecular decay or infrared 
emission. 

3.1. Vibrational redistribution 

Neglecting the influence of the excited states, collisions and 
UV-visible photon absorption under astrophysical conditions 
leads to an increase of the internal energy of the cluster, which 
can then be dissipated either by evaporation and/or IR emis- 
sion. The energy redistribution between inter- and intramolec- 
ular vibrational modes is required for the cluster to be con- 
sidered in thermal equilibrium upon nucleation of a molecule. 
Similarly, the absorption of a photon by a single molecule in 
the cluster will tend to heat the intramolecular modes of this 
specific molecule first by internal conversion, before the ex- 
cess energy flows into intermolecular modes, and subsequently 
into the intramolecular modes of the other molecules. The time 
scales involved in these fundamental redistribution processes 
should be considered as reference data in the next stages of our 
investigation. 

As a first application of our model for PAH clusters, we 
have simulated the energy redistribution of intramolecular en- 
ergy toward intermolecular modes by heating a single molecule 
in the cluster and monitoring the amount of kinetic energy in all 
other molecules as a function of time. The kinetic energy pro- 
vides a direct estimate of the vibrational or thermal energy. We 



consider redistribution in the (coronene)5 cluster as a typical 
example. This cluster consists of a short stack, for which either 
the outermost or the central molecule are chosen to carry the 
initial thermal excitation. The starting configuration of the clus- 
ter is taken as the lowest-energy structure, and no extra kinetic 
energy is given to the intermolecular modes. To quantify the re- 
distribution of energy into vibrational modes of all molecules 
in the cluster, we introduce the root mean square fluctuation Sj 
of the kinetic temperatures of the molecules in the cluster: 



^T{t)^-^.y\{Ti{t)-{T{t))f 



(2) 



where r,(f) is the average kinetic temperature restricted to 
molecule / and {T(t)) is the usual average kinetic temperature 
of the entire system. 5r(0 provides an instantaneous measure 
of the standard deviation of the individual temperatures of all 
molecules in the cluster Initially this quantity is large, since 
only one molecule carries kinetic energy. Eventually the clus- 
ter reaches equilibrium and the temperatures of all molecules 
in the cluster are identical. The rate at which dTif) decreases 
yields an estimate for the time scale r of energy redistribution 
in the system, according to ^^(f) ~ 6T{Q)e^'^^- 

Figure|4]shows the variations of 6t as a function of time for 
(C24Hi2)5 in which the initial excitation is localized on differ- 
ent molecules. The decrease of 6t is indeed exponential, with 
decay time constants approximately equal to r = 1.47 + 0.09 
ps when the central molecule is vibrationally excited, and t = 
2.67 + 0.28 ps when the outermost molecule is excited. These 
time scales are of the same magnitude as the softest vibrational 
modes of isolated coronene molecules (>Martia.l996.) . This is 
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Fig. 4. Root mean squares 6t of the kinetic energies of the 
molecules in a cluster of five molecules. At f = 0, the initially 
kinetic energy excess (13 eV) is localized in the outermost (cir- 
cles) or central (squares) molecule within the stack. The dashed 
lines are linear fits. 



not surprising, because the slowest modes constitute a lower 
limit to the rate at which energy is redistributed homogeneously 
within all internal degrees of freedom. Finally, the faster redis- 
tribution when the central molecule is heated is expected due 
to the presence of fewer intermolecular modes involved. 

3.2. Evaporation 

Coronene clusters could a priori dissociate into entire 
molecules, or via intramolecular fragmentation. The loss of a 
molecu le from a cororiene clu ster [typically 1 eV of binding 
energy feapa cioh et al. P2Q0.5a')1 is much more likely to happen 
than t he breaking of chem ical bonds, which requires more than 
4 eV j-TobUn et alJl2Q06l). In agreemen t with the experimen- 
tal conclusions ofl Schmidt et al.l I 2006"), we thus only consider 
here the dissociation corresponding to single molecule evapo- 
ration. 

Evaporation rates of molecular clusters are not available by 
direct molecular dynamics simulations except at very high ex- 
citation energy. When the extra energy deposited in the cluster 
is marginally higher than the dissociation energy, fragmenta- 
tion is a rare event, which can take place over time scales much 
longer than the typical times that are reachable by conventional 
molecular simulations. Statistical rate theories are much more 
convenient for calci ilating such rates. The phase space theory 
(PST) developed bv lChesnavich & BowersI lll977h provides an 
accurate framework for unimolecular dissociation. In previous 
work, two of us (.Calvo & Parneix 2004) showed that PST pre- 
dictions were in quantitative agreement with numerical simu- 
lations for molecular systems, provided that the rotational and 
vibrational densities of states (DOS) were correctly evaluated. 



3.2.1 . Phase space theory 

The success of PST is essentially due to the rigorous conserva- 
tion of energy and angular momentum, and to the inclusion of 
anharmonic vibrational effects. Here we are mainly interested 
in the determination of absolute rate constants for the evap- 
oration of a single molecule fro m a thermaliz ed cluster PST 
expressions for the rate constant ll.Tarroldll99li) contain a set of 
factors that are not always obvious to estimate with accuracy. 
Instead of assuming specific values for these unknown factors, 
we follow Weerasinghe & Am^ ([1993) and use the results of 
molecular dynamics simulations at high internal energy to cal- 
ibrate these unknown factors. 

Performing molecular dynamics simulation of the cluster 
with both inter- and intramolecular interactions is not realis- 
tic in the context of thousands of trajectories over time scales 
governed by the intermolecular motion. Therefore we had to 
simplify the modeling by freezing intramolecular motion. With 
this assumption, the molecules are treated as rigid bodies, and 
larger time steps and algorithmic improvements such as the use 
of quaternion coordinates allow an efficient simulation of the 
evaporation process. Even though the rigid body approxima- 
tion may not be physically relevant at high energies, it cali- 
brates the statistical theory, that we will subsequently use in a 
lower excitation regime. 

In this subsection, energies are thus measured on the inter- 
molecular potential energy surface //inter only. It is not our goal 
here to describe the PST form afism in great details, and th e 
reader is refered to the work bvlCh^ ^navich & BowersI i ll 977h . 
IWeerasinghe & Amad (Il993h . ICalvo & Parneix (.2004) for fur- 
ther information. The difi'erential rate of unimolecular evapora- 
tion from a molecular cluster Xn+\ X„+X (here X=C24Hi2), 
at total rovibrational energy E and angular momentum J, and 
leading to the loss of kinetic energy Str into translation St and 
rotation is expressed as 



R(E,J,etr) = Co 



Q.'(E-Eo-s,,)r{J,s,,) 

□(£■ - Er) 



(3) 



In the above equation Eq is the dissociation energy, that is the 
binding energy difference between the parent and the product 
clusters, E^ is the parent rotational energy, equal to BJ- where 
B is the rotational constant. Q and Q' are the vibrational den- 
sities of states of the parent and product clusters, respectively. 
F' is the rotational density of states of the two products, ob- 
tained from the mechanical constraints on angular momentum 
at a given translational-i-rotational energy. An estimation of F' 
requires treating the products as rigid bodies. Co is a constant 
that we avoid calculating by performing high energy molecular 
dynamics simulations for calibration. The rate constant k(E, J) 
is obtained from the dififerential rate R(E, J, en) by integration 
over the available range of St, . 

The ingredients of the PST calculations are threefold. First, 
the vibrational densities of states Q and Q.' are calculated us- 
ing parallel tempering Monte Carlo simulations ('Neir otti et alJ 
.2000) analysed with a multiple histogram method. These calcu- 
lations consist of running independent Monte Carlo trajectories 
or replicas at different temperatures, occasionnally attempting 
to swap configurations between adjacent replicas. Here we use 
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50 different temperatures in the range 10 K-1500 K distributed 
according to a geometric progression. For each trajectory, 6 10^ 
Monte Carlo cycles (1 MC cycle =« individual steps) are car- 
ried out for the calculation, the first 10^ steps being discarded 
for equilibration. 

The rotational density of states F' requires integrating the 
number of rotational states x with specified internal angular 
momentum of the products J, - J\ + J2 and specified orbital 
momentum L under the constraint that the translational energy 
released e, is larger that the centrifugal barrier e' . The calcula- 
tion becomes significantly simpler by assuming that the initial 
angular momentum is small or, more precisely, that the rota- 
tional energy is small with respect to the vibrational energy. 
In that case Jr = L in modulus, and the integration is one- 
dimensional over the range of available values of Jr'. 

r(e„.,y-o)= f' x\s:,J,-)dJr- (4) 
Jo 

In the previous equation e* = etr - with eHjr) the energy of 
the centrifugal barrier 

The function y to be integ r ated h as been calculated ex- 
actly by IChesnavich & BowersI ( 1 19771) in a number of situa- 
tions for products with various shapes. However, coronene is 
a highly oblate molecule, and expressions are only available 
if the other product has a high symmetry. Therefore we had 
to restrict our study to sizes for which the product cluster is 
roug hly spherical. From ou r previous results on cluster struc- 
ture jRapacioli et alJEoOSah . two sizes fulfill this criterion to 
a good approximation, namely n - 3 (B » 3.8 10"^ cm"') 
and n = 12 (B » 1.4 10"* cm"'). Thus our investigation of 
unimolecular dissociation will focus on the two parent clusters 
containing 4 or 13 molecules. 

As a third ingredient, the PST calculation also needs an ex- 
pression for the dissociation potential between the products, in 
order to locate the centrifugal barrier and estimate its height 
e ' . At long range r — > 00 the interaction between a cluster of 
PAH molecules and a single PAH is essentially of the disper- 
sion form, which we correct at short distances using the ex- 
pression -Ce/{r - ro)^. The values of Ce and ro were obtained 
from a series of Monte Carlo simulations with the distance be- 
tween the two products constrained to increasing values. We 
find Ce = -4.3617 10* kJ A^/mol (resp. Q = -1.6574 10^ 
kJ A*/mol) and ro = 2.27 A (resp. ro = 3.73 A) for n = 4 (resp. 
n = 13). 

3.2.2. Calibration 

The statistical theory described previously can be used to calcu- 
late kinetic energy release (KER) distributions, given by Eq. ^ 
after normalization over Str- The KER has the advantage over 
the rate constant that it does not include the unknown factor 
Co, hence it provides a good testing ground for the accuracy of 
PST. 

For each of the two parent clusters, we performed a sam- 
ple of 10 000 molecular dynamics simulations at relatively high 
excess energy. For each trajectory, the cluster was initially pre- 
pared cold (10 K), and was suddenly heated. The simulation 



was ended when dissociation of a single molecule occured 
within 100 ps. Once evaporation was detected, the dissociation 
time was taken as the latest moment when the radial velocity 
of the dissociating molecule was pointing inward. Here the MD 
results can be considered as numerically exact, and as a refer- 
ence for assessing the PST approach. 

The distributions of KER obtained from the MD simula- 
tions are compared with PST predictions in Fig.|5]for the dis- 
sociation of (C24Hi2)4 at 2 eV excess energy. 
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Fig. 5. Probability distribution of the kinetic energy released 
during the unimolecular evaporation of (C24Hi2)4 at 2 eV in- 
termolecular energy obtained from MD simulations and from 
PST based on several approximations (see the text for details 
about these approximations). 

The PST calculations were performed at three different lev- 
els of approximation. In its most rigorous form, the vibrational 
density of states is the anharmonic form extracted from Monte 
Carlo simulations, and the rotational density results from in- 
tegrating Eq. (|4j above. The harmonic approximation for the 
vibrational density is simply Q.(E) oc E^"' with v = 6n - 6 
the number of independent degrees of freedom in the system (3 
rotational and 3 translational modes per molecul e). The Klots 
approximation at low angular momentum [see llKIotslll97lt 
ICalvo & Parneixll2004 for more details] is similarly r'(eti) °^ 
ej^ with a = (r - l)/2, r = 6 being the number of rotational 
degrees of freedom in the products. 

From Fig.|5]we generally find a very good agreement be- 
tween the MD results and the PST predictions. In the case of the 
tetramer, the harmonic approximation for the vibrational DOS 
performs well. This is consistent with the existence of only a 
single (stacked) isomer for this cluster at the energies consid- 
ered here. The Klots approximation for the rotational density 
does not affect the results quantitatively. This indicates that the 
results should be poorly sensitive to the exact shape of the prod- 
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ucts. For the larger cluster (distributions not plotted) the agree- 
ment remains very good but anharmonicities are found to be 
more significant. 

The present agreement between molecular dynamics data 
and the PST results is encouraging and allows us to fit the miss- 
ing parameter Co to reproduce the rate constant obtained from 
MD. Subsequently the statistical properties can be calculated 
from PST for any excitation energy. This is especially interest- 
ing at low energies, for which molecular dynamics simulations 
are not practical. 
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-0.2084 
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2.4862 


-0.1666 


0.00449 



Table 1. Parameters of the polynomial fits of Etot, as a function 
of iiintra per molecule (ak) and Winter (bk) for the (coronene)4 
and (coronene)i3 clusters. 



3.2.3. Evaporation rates 

The variations of the rate constant with internal energy are rep- 
resented in Fig.0for the two parent clusters. The physical con- 
tents of these results will be discussed later in the light of other 
competing decay channels. 
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4. Competition between evaporation and IR 
cooling 

Infrared emission is a possible relaxation energy channel for 
an excited cluster and can be in competition with evaporation. 
While evaporation mainly involves the intermolecular energy, 
most of the IR emission is driven by intramolecular modes. 
Therefore, in order to compare the evaporation and IR emis- 
sion rates, the inter- and intramolecular energies must be first 
related to each other 



4.1. Energy distribution 

In the next subsection, the spontaneous IR emission rates are 
calculated for each molecule against their intramolecular en- 
ergy /iintia- The evaporation rates as a function of intermolec- 
ular energy isinter are given in Sec. 13.2.51 fiintra and /Simer both 
have to be expressed as a function of the total cluster energy 
Zitot. We first assume that the PAH cluster remains at thermal 
equilibrium. Because most intramolecular mode s have a high 
frequency, the zero-point effects are important tiSchmidt et alJ 
2006). The intermolecular modes, on the other hand, are all 
essentially soft and the intermolecular energy is equally dis- 
tributed among those modes following classical equipartition. 

We calculate Eintm explicitly as a function of tempera- 
ture from the known vibrational frequencies {w,) of the single 
coronene molecule: 



EintmiT) = ncUi + 



1 



exp{Tia>i/kET) - 1 



(5) 



Fig. 6. Intramolecular energy per molecule and total inter- 
molecular energy as a function of the total energy in the cluster 
for (C24Hi2)4 (top panel) and (C24Hi2)i3 (bottom panel), after 
removal of harmonic zero-point energies. 

An independent validation of phase space theory is ob- 
tained by comparing the rates from sets of molecular dynam- 
ics simulations performed at different energies to the PST pre- 
dictions. We thus performed additional simulations for inter- 
molecular energies in the 1.5-1.9 eV range for (C24Hi2)4 and 
at 5.5 eV for (C24Hi2)i3. The results shown in Fig.0are in very 
good agreement with the statistical theory, even though a log- 
arithmic scale has been used for the graph. The discrepancy is 
found to be at most a factor 1.85, which is satisfactory for a rate 
constant known to span orders of magnitude. 



being the Boltzmann constant. The frequencies {w,) are ob- 
tained from the DFT calculations of Martin ( 199.6) (which are 
more accurate than the present TB model). 

Similarly, we calculate isintei as a function of tempera- 
ture using the intermolecular frequencies from R apacioli et alJ 
(2005a). The total energy of the cluster is the sum of the intra- 
and intermolecular energies. Figure |6lrepresents the variations 
of the intramolecular energy per molecule fiintra and of the to- 
tal intermolecular energy isinter against Etot- These curves were 
fitted using polynomial expressions of the form: 



Etot - ^ flj(:(£'intra)* 
= 2/7,(£i„,ei)' 

with fitting coefficients given in Table [fl 



(6) 
(7) 
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Fig. 7. Evaporation rates obtained with PST (solid lines) and MD simulations (symbols) as a function of the internal energy for 
coronene clusters with 4 (left panel) and 13 (right panel) molecules, respectively. The IR emission rates (Joblin et al. 2002) are 
also represented as dashed lines. The total, intra- and intermolecular energy scales are related to each other through Eqn. (|6j and 
0. 



4.2. Spontaneous IR emission rates and evaporation 
rates 

We address here the estimation of IR emission rates for neutral 
coronene clusters. In a first approximation, the effects of clus- 
tering on the IR emission of a specific molecule in the cluster 
will be neglected. Thus the IR emission rate of the full clus- 
ter is obtained, considering the rates of individual molecules. 
The latter quantity explicitly depends on the energy /Sintra- 
Calculations of the IR emission of coronene molecules in in- 
terst ellar conditions we re reported by Joblin et al. (2002) [see 
also jMulas e t al. 2006) for a more refined version of this ap- 
proach]. In this Monte Carlo code, the IR emission rates are 
described with microcanonical statis tics using an exact count- 
ing method for the density of states dBever & Swineharll973l) 
based on harmonic oscillators. 

To discuss qualitatively the competition between the evap- 
oration and IR emission processes, a first convenient step is 
to define a mean emission rate. This is achieved by summing 
the emission rates of all IR active vibrational modes of all 
molecules in the cluster : 

^IR(£,„tra) = Z Z ^IR (^y)(^."tra), (8) 
n i,j 

where k^^^~^(vj) is the spontaneous emission rate for a vj IR 
mode in the v, ^ v, - 1 transition at the available energy £intra- 
The variations of the mean IR emission rates ^ir and the 
evaporation rates ^evap calculated in section ll!231 for (C24Hi2)4 



and (C24Hi2)i3 are represented in Fig.0against fitot- The clus- 
ters are assumed to be at thermal equilibrium, the energy repar- 
tition of section RTTI was used. 

A critical energy E* can be defined as the energy of equal 
rates, A;evap(£*) = ^ir(£*)- This energy is 10.4 eV for (C24Hi2)4 
and 18.6 eV for (C24Hi2)i3. 

While the evaporation rates increase very strongly (expo- 
nentially) with Etot, the IR emission rates show much smaller, 
polynomial-type variations. This suggests that the probability 
of evaporating a cluster before IR cooling strongly decreases 
for energies smaller than E* and becomes maximum for ener- 
gies larger than this threshold. 

5. Modelling the photophysics of a neutral 
coronene cluster in a radiation field 

5.1. Coronene cluster in a radiation field 

The Monte Carlo calculation of 'joblin et al.' ('2002") is based 
on the exact s tochastic method given by Gillesp ie ( 1978) and 
lBarkei< ill 9831) . The method consists of propagating 'trajecto- 
ries' that are Markovian random walks over the energy lev- 
els of the molecule. These trajectories depend on the follow- 
ing events: absorption of UV photons, e mission of IR photon s, 
but also fragmentation, as described in LlobUn et alJ ('2006') . 
We have adapted this model to describe the photophysics of 
coronene clusters in a specific radiation field. The absorption 
rates were calculated using the absorption cross section cr of 
the coronene cluster and a description of the astronomical radi- 
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ation field. Since cr is not known, we have scaled the cross sec- 
tion of the coronene molecule measured by Joblin et al. ( 19921) 
by the number of molecules. The IR emission rates as a func- 
tion of iiintra wcrc described in Sec. 14.21 The evaporation rates 
were estimated by fitting the curves of Fig. Q with Winter calcu- 
lated from Etot, as explained in section BTTI 

The model allows us to calculate the time evolution of the 
internal energy and the associated events (photon absorption, 
emission, or fragmentation). Simulations are run over a given 
number of absorbed photons (10^ and 10'' for the lowest and 
highest absorption rates, respectively). 

5.2. Simulation of NGC 7023. The muitiphoton 
absorption phenomenon 

The photophysical evolution of (C24Hi2)4 and (C24Hi2)i3 clus- 
ters has been simulated in the environment of the northern 
photodissociation region (PDR) of NGC 7023. According to 
[Rapacioli et al. ( 2005b), different regions corresponding to var- 
ious depths inside the molecular cloud were selected. The depth 
is quantified by the selective extinction A i/. 

The radiation field in the different r egions was cal- 
culated using the Meudon PDR code llLe Bourlot et alJ 
( 1 1993b . iLe Petit et all t.OO&i\. The inputs for this calcula- 
tio n were modifi e d with respect to previous calculations 
of iRapacioli et a l] il2005bl) . For the total-to-selective extinc- 
tion, we used a value for Ry of 5.5 according to the re- 
cent study by WilLgt al, (2006). The stellar spectrum used 
for HD200775 is from the Ku rucz library with a tempera- 
ture of 1 5 000 K ( lKurucz"l99l'). The higher temperature de- 
rived by Ivan den Ancker et al. (1997) was found to be in- 
consistent with th e VUV flux o bserve d on the star by the 
lUE satellite [see 'Valenti et al.' ('200Cf) and lUE archive at 
http://www.vilspa.esa.es/iue/iue.htmll. The radius of the star 
was taken as \QRQ which is typi cal for a Be star . We performed 
the same analysis as in Rapaciol i et alJ ( l2005bl) and inferred a 
hydrogen density of 7 10^ cm""* and a value Ay - 0.41 for the 
region of maximum emission in the H2 9.7 /im band. 

The radiation field was calculated for Ay values of 0.41, 
1.0, 1.6, 2.5, 3.5, 5.0 and 6.4. The variations of the correspond- 
ing UV absorption rates are represented in Fig.|8] 

The average timescales of evaporation are shown in Fig. 
|9]for (C24Hi2)4 and (C24Hi2)i3. While the small evaporation 
timescale for (C24Hi2)4 can be explained by the low critical 
energy E* of this cluster (see section 14. 2> . the evaporation 
of the larger (C24Hi2)i3 has its critical energy above the cut- 
off of the radiation field. This is illustrated in Fig. [TO] which 
displays the distributions of internal energy of the clusters 
prior to their evaporation. The mean energy for dissociation 
is in the range 10.5-11.5 eV for (C24Hi2)4 and -14.5 eV for 
(C24Hi2)i3. This value is slightly shifted towards lower ener- 
gies at Av = 5.0 because of the extinction, which preferentially 
affects the more energetic photons. Another effect that can be 
observed in Fig. ^| is the absence of the high-energy tail for 
Ay - 5.0. Such events are only possible if the molecule can ab- 
sorb two VUV photons without totally cooling by IR emission 




Fig. 8. Absorption rates of C24H12 at several depths inside the 
northern PDR of NGC 7023 : Ay^ 0.41, 1.0, 1.6, 2.5, 3.5, 5.0 
and 6.4 from top to bottom. For clusters, the absorption rates 
should be multiplied by the number of molecular units. 
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Fig. 9. Characteristic timescale for the evaporation of 
(C24Hi2)4 (circles) and (C24Hi2)i3 (squares) as a function of 
Ay in the northern PDR of NGC 7023. 

between the two absorptions, which requires in turn requires 
sufficiently high absorption rates. 

From Fig. |9] it is clear that (C24Hi2)i3 clusters have ac- 
cumulated the energy of more than one photon to reach the 
dissociation threshold, even in the more embedded regions. 
Figurefm shows the internal energy contained in the clusters 
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Fig. 10. Normalized distributions of the internal energy in 
(C24Hi2)4 and (C24Hi2)i3 clusters prior to their evaporation 
(black and white bars respectively). The calculations were per- 
formed in two different radiation fields at {a) Ay = 0.41 and (b) 
Av = 5.0. 




a. 10 
i 10-3 

LU 

10-^ 



Q 1 23456789 1011 1213 2 4 e 8 3 12 14 16 18 

E„ (eV) E (eV) 



Fig. 11. Internal energy distributions of the clusters before 
absorption of the VUV photon leading to dissociation of 
(C24Hi2)4 (left) and (C24Hi2)i3 (right). The calculations were 
performed in two different radiation fields at (a and b) Av = 
0.41 and (c and d) Ay - 5.0. The distributions are normalized 
to the total number of absorbed UV photons. 



before absorption of the VUV photon leading to dissociation. 
In the case of the large clusters (13 units), a clear threshold is 
seen at 3 eV. This threshold is very sharp at low energies, as 
could be anticipated: there is no event in the 2.85 eV channel 
of the histogram and 60-70% of the total number of events are 
found in the 3.0 eV channel. This fraction drops very quickly 
at higher energies with ~30% in the 3.25 eV channel and the 
remaining 1-10% spread over channels up to 15.5 eV and 
13.5 eV for Al/ of 0.41 and 5.0, respectively. 

The peak in the 3.0 eV channel can be interpreted as the 
consequence of truncating the mean energy stored in the clus- 
ters at the threshold where the clusters are in equilibrium with 
the radiation field. Indeed, the model shows that a full dissipa- 
tion of the energy absorbed by IR emission cannot be achieved 
before another photon is absorbed. At low internal energies, the 
cooling efliciency decreases as only low frequency modes with 
weak intensities are involved. 

Contrary to large clusters, (C24Hi2)4 species can be evap- 
orated by a single VUV photon, multi-photon events playing 
only a minor role, as seen in Fig.^2 In this case, the eV 
channel corresponds to clusters that evaporate as soon as they 
absorb their first UV photon. The peak around 0.8 eV corre- 
sponds to the mean energy stored in equilibrium with the radi- 
ation field. In this case, we can conclude from the observation 
of a population in the eV energy channel that accumulation of 
energy is not required for the evaporation of (C24Hi2)4 clusters. 
This is obvious since energies lower than 13.6 eV are already 
required. 

6. On the competition between formation and 
destruction of PAIH clusters under interstellar 
conditions 

6.1. Growth constants 

We estimate here the rate constant of the nucleation process 
associated with the growth of a single molecule on a cluster of 
/ - 1 molecules, both species being neutral. In the interstellar 
medium, the molecular gas can be considered as very rarefied, 
hence molecules and clusters have enough time to release their 
internal energy via IR emission before any collision occurs. In 
addition, we will assume that the molecular gas is thermalized, 
Boltzmann laws governing the distribution of translational and 
rotational energies of both clusters and single molecules. 

If dn, represents the number of clusters of size / formed 
upon the nucleation reaction 1 -H (/ - 1) — > / per unit of time dt 
and inside a unit of volume, the growth constant ki is defined 
such that: 



Shi 
St 



(9) 



In this equation «! and n,_i are the densities of single molecules 
and clusters of size /- 1 per volume unit, respectively, kj can be 
expressed from standard collision theory cLandau et al_ 198ll) 



as: 



ki^ JJJ /i(vi)/;-i(v,-,)|K 



Xpi,^(yr;Er)d\ld\i-lfiEr)dEr. 



(10) 
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In the above equation, /i(vi) and ./;-i(v,_i) are the probability 
densities of the incident molecule and the cluster, which de- 
pend on their respective velocities Vj and v,_i only. f{E^) is the 
probability distribution of the rotational energy available at the 
beginning of the collision, ctq is the nucleation cross section, 
taken here as the geometrical cross section. Vr = Vi - Vi_i is the 
relative velocity vector. pi+(Vr',Ey) is the probability that the 
collision leads to the formation of cluster size /. Assuming that 
/i(vi) and //-i(v,_i) both follow a Boltzmann law at tempera- 
ture T, Eq. ilQ\ can be written more explicitly as: 



Similar calculations have also been performed with the proba- 
bility /:>, _ of destroying the parent cluster The corresponding 
rates are found to be a few orders of magnitude smaller than 
the growth rates, hence we conclude that this process can be 
neglected. 



X exp 



J J l|VrllP/,+ (Vr; 



(OTi H-OT,_i)v^ +yL<V? 



E,)a-Q 



(11) 



In the latter equation, nii and m,_i are the masses of the inci- 
dent molecule and the cluster, respectively; Vc is the velocity 
of the centre of mass of the entire system. The integration of 
this equation over Vc gives: 



(TQPiAEu E,)e-^^"'^^dEJiE,)dEr, (12) 



where E^ is the translational energy in the centre of mass refer- 
ence frame, and // the reduced mass of the system. 

In order to calculate quantitative values of k,, the collision 
cross section is approximated as the geometric area of the PAH 
molecule (7.1 10"'^ A"^). In section fT/Jl we only calculated 
values of pi^+iE^; E^) for the specific case £, = E^. However a 
lower bound to ^, can be estimated as follows. 

At a given translational energy Ei, the cluster growth prob- 
ability decreases with the increase of rotational energy E^, as it 
implies an increase in the total energy. Hence 



PiAEuE,>E^)<PiAEuE°). 



(13) 



Similarly, at a given Ei-, the cluster growth probability de- 
creases with increasing E^: 



PiAEt>E°;E,)<pi,^iE°;E,). 



(14) 



These two trends have been checked using explicit atomic sim- 
ulations, and performing a few tens of molecular dynamics tra- 
jectories of the nucleation process. A lower bound to ki can be 
calculated assuming that 



PiAEt; E,) 



■piAEuE,) if£, >£,.: 
.PiAEx;E,) ifEi<E,. 



(15) 



Moreover, the collision constants /:>, + provide an upper bound 
to the growth constants, by taking pj + = 1 in Eq. il2i . 

Figure shows the lower and upper limits to kj with 
i = 2,3,4 obtained in the temperature range 10^-10"* K. The 
two limits appear to be very close. This is due to the fact that, 
except for extremely high coUisional energies, nearly all col- 
lisions lead to effective growth. In a first approximation, the 
growth rate can then be taken as the collision rate. We also note 
that the growth constant increases monotonically with temper- 
ature. Indeed a temperature increase enhances the rate of colli- 
sions which, for T < 2000 K, nearly always lead to nucleation. 
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Fig. 12. Collision rates (dashed lines) and growth rates (solid 
lines) for the nucleation processes l-i-l— >2, 1h-2— »3 and 
1 -H 3 — > 4 plotted as a function of temperature. 



6.2. Competition between formation and destruction in 
neutral PAH clusters 

In their study of the reflection nebula NGC 7023, 
iRapaciol i et al. (2005b) gave hints that PAH molecules 
are likely to be condensed as clusters. These clusters are 
evaporated at the border of the cloud where isolated PAHs 
are observed. A minimum size of 400 carbon atoms for the 
observed clusters was also inferred from the observations. 

An order of magnitude for the typical timescale of PAH 
clustering in this reflection nebula can be estimated from the 
calculations reported in section |6l We use a density riu - 
7 10^* cm"-' for NGC 7023, as discussed previously (sec. l5.2>. 
and a carbon abundance ric - 2 lO^^ne following Sofia et alJ 
( I2OO4I) . We consider that about 20% of the carbon is contained 
in PAH molecules. For simplicity, and in order to apply our 
previous results, we also assume that all PAH molecules are 
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coronene. This leads to a density of coronene molecules of 
1.12 10"^ cm"^. According to the PDR model, the kinetic tem- 
perature varies from 500 K to 60 K in the range - 0.41-6.4. 

We find that the aggregation timescale is typically lO'*- 
8 10"* years. This is much longer than the photoevaporation 
timescales derived for (C24Hi2)4 and (C24Hi2)i3 in section ls!2l 
which are below 1.5 10^ hours (17 years). Therefore, these 
clusters should be destroyed by the UV flux much faster than 
they can be reformed. This also implies that, when a large clus- 
ter starts to evaporate, the small fragments are also likely to 
photoevaporate themselves. It can be deduced from this theo- 
retical work that, in the case of NGC 7023, the minimum size 
of the coronene clusters should be larger than 13 molecular 
units (312 carbon atoms), which is consistent with the mini- 
mum size of 400 carbo n atoms inferred from the observations 
jRapacioli et al.l2005bl) . 

6.3. Effects of ionization 

Until now we have focused on the description of neutral clus- 
ters. Our model can be used to constrain the size distribution of 
PAH clusters as a function of depth in the molecular cloud. One 
would like to be able to constrai n the charge as ^yell, similarly 
to wh at is done for sing le PAHs ( Bake s et alJ200lllLi & DraTiij 
I2OO2I) . While we can give hints about charge effects, funda- 
mental studies are clearly required before considering further 
modeling. 

Even though PAHs are likely to carry a charge in 
many objects of the ISM (iHudgins & A llamandola 199^ 
lAUamandola et al.' 'l989'). R apacioli et aTn2005bft have shown 
that the different populations of very small dust particles are 
spatially decorrelated in two PDRs NGC 7023 and pOph-SR3. 
As the distance from the stars increases, the mid-lR emission is 
found to be dominated first by positively ionized PAHs, then by 
neutral PAHs, and finally by VSGs. Similarly to single PAHs, 
PAH clusters could carry a charge. We expect this charge to be 
most likely positive, since the lifetime of anions will be limited 
by the easy photodetachment in these UV-irradiated regions 
where PAH clusters can be observed by their mid-lR emission. 
Charge effects can affect both the formation rate and the evap- 
oration rate for different but related reasons. 

Let us first consider the case of multiply charged systems. 
The nucleation between two charged subsystems seems very 
unlikely due to the strong Coulomb repulsion. We thus ex- 
clude the reactions of positively charged clusters with PAH 
cations. The stability of an existing assembly of PAH molecules 
will critically depend on its charge state, in addition to its 
size. Other factors, including the size of the PAH itself and 
the excitation mechanisms, presumably play an important role 
as well. In a multiply ionized cluster. Coulomb fragmenta- 
tion will generally enhance the possibility of emitting sin- 
gle (charged) molecules. In addition, new fragmentation chan- 
nels will be open if a single PAH molecule carries more than 
one charge, and various intermolecular and/or intramolecu- 
lar Coulomb multifragmentation processes could take place. 
While the photofragmentation dynamics of a PAH assembly 



(or even a single PAH cluster) would be worth investigating in 
its own, it is beyond the scope of the present paper. 

Based on previous theoretical calculations on small singly 
charged clusters of aromatic molecules such as benzene, naph- 
talene or anthracene (Piuzzi et al. 200 2: Bouvier et al. 2 002t) . 
cationic clusters are expected to be significantly more sta- 
ble than their neutral counterparts from the energetic point of 
view. This increased stability reflects the charge delocalization 
over several molecules in the cluster, a collective effect which 
cannot be described using simple explicit potentials. Because 
the evaporation rate constant varies approximately exponen- 
tially with the dissociation energy, cationic clusters will be sig- 
nificantly more resistant to the decay of a neutral molecule. 

In the region where PAH clusters are observed, P AHs are 
expected to be mainly neutral llRapacioh et alJ2005bl) . and we 
can limit our discussion to the case where a neutral molecule 
collides with a (PAH),^ cluster. For the nucleation rate, the PAH 
molecule being nonpolar, the long-range ion/neutral interac- 
tion scales as -C4/r* where the constant C4 is proportional 
to the electronic dipole polarizability a of the neutral PAH. 
In a very simple approach, we use a Langevin model to esti- 
mate the collision rate constant. Based on our previous results 
on neutral/neutral collisions, we can reasonably assume that all 
collisions between a neutral PAH and a cationic cluster lead 
to an effective growth of the cluster The Langevin prediction 
for the nucleation rate is t hen g iven by k„ur] - 2neq{alij)^^^ 
[see for instance iHerbstI (l200lh l. where e is the electronic 
charge, q the ionization degree iq - \ for singly-charged clus- 
ters) and the reduced mass of the reactants. An average po- 
larizability of the coronene molecule, as estimated from the 
components of the polarizability tensor via density-functional 
theory calculations (Malloci 2006), is a = 46 lO"^"* cm-'. 
Using this value the nucleation rate is approximately A;nuci — 
9 10 '°/^ cm^s Together with the PAH density derived pre- 
viously for NGC 7023, this leads to aggregation timescales of 
[2.0-3.5]/^ 10^ years. 

For the specific case of NGC 7023, considered in section 
16.21 the conclusion that small clusters are destroyed faster than 
they can be reformed could thus be modified if the clusters 
carry one single positive charge, as this would decrease their 
evaporation rate (mostly through increasing their binding en- 
ergy) but also enhance the nucleation rates. 

6.4. Astrophysicai implications 

As seen above, a single positive ionization is expected to sta- 
bilize PAH clusters with respect to evaporation, and to favor 
their formation by nucleation of a neutral molecule. Our con- 
clusion that neutral PAH clusters are evaporated efficiently in 
the northern PDR of NGC 7023 without being reformed could 
thus be revised in the case of cationic clusters. 

The linear scaling of the growth constant with inverse den- 
sity suggests another mechanism that could affect the condi- 
tions of reformation. The northern PDR of NGC 7023 is hetero- 
geneous, with some high-density filaments (n ~ 10^-10* cm"^) 
embedded in a more diffuse medium with n ~ 10"* cm"-* 
LFuente et al. C2000i) and references therein]. With a value of 
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5 10'' cm"'' for the density, the timescale for formation could 
be lowered to 140-700 years for a neutral cluster, and to SO- 
SO years for an ionized cluster Considering the dissociation 
rates reported in Fig.|5]for neutral species, and taking into ac- 
count the extra stability of the cationic species, the PAH clus- 
ters could thus survive much longer and closer to the border 
of the cloud. A scenario in which the growth of PAH clusters 
occurs in the PDR itself is however unlikely because it would 
involve the survival of dimers in PAH-emitting regions. This 
suggests that such clusters would be formed in the UV-shielded 
parts of the cloud. There is time to form them during the cycle 
of matter from diffuse to dense clouds. Indeed the quiet phase 
of molecular clouds will typically last 10^ years, allowing clus- 
ters of tens of molecules to be formed. The formation of young 
stars in these clouds creates PDRs, which are fed in with clus- 
ters and rapidly submitted to photoevaporation. Dynamical ef- 
fects in PDRs such as an advancing photodi ssociation region 
(lLemaireetal .'l999: Fuente et al. 1999,'2000|) could then bring 
these clusters into the UV-irradiated regions. 

7. Conclusion 

The formation and physico-chemical evolutio n of PAHs is a 
ques tion of great interest in astrochemistry. Rapaciol i et alJ 
ll2005b.) presented strong evidence that PAH clusters are good 
candidates for small carbonaceous grains. These authors sug- 
gested that those grains are indeed PAH clusters, and a min- 
imum size of 400 carbon atoms per cluster was inferred in 
NGC 7023. In this paper, the competition between the forma- 
tion and destruction processes of such PAH clusters in the neu- 
tral state was investigated theoretically under realistic astro- 
physical conditions. Coronene clusters have been considered 
here as representative prototypes of PAH clusters. 

The growth of neutral PAH clusters by nucleation of indi- 
vidual molecules was simulated using a model combining an 
explicit atomic force field to describe the intermolecular inter- 
actions, along with a quantum tight-binding approach for the 
intramolecular interactions. It was found that, under interstellar 
conditions, most of the collisions lead to cluster growth. This 
was interpreted as due to the large number of intramolecular 
modes that contribute to absorbing the excess collision energy. 

Evaporation was treated in the framework of phase space 
theory, supplemented with molecular dynamics simulations. 
The clusters (C24Hi2)4 and (C24Hi2)i3 were investigated, and 
a good agreement between the simulation results and the pre- 
dictions of the statistical theory was found for the distribu- 
tion of kinetic energy released after evaporation. The calculated 
evaporation rates were then introduced in an infrared emission 
model based on microcanonical statistics. This model allowed 
the simulation of these clusters in a specified interstellar radia- 
tion field. 

When applied to the reflection nebula NGC 7023, our 
model indicates that (C24Hi2)4 and (C24Hi2)i3 are both photo- 
evaporated much faster than they can be reformed. In the case 
of (C24Hi2)i3, this evaporation is possible if the UV absorption 
rate is faster than the IR emission rate, so that at least ~3 eV 
of internal energy are stored in the cluster before absorption 
of a VUV photon. These results provide a plausible explana- 



tion for the absence of small VSGs observed in this PDR. They 
also suggest that PAH clusters are formed inside the molecular 
cloud and transported in the PDR by dynamical effects. 

The present investigation is a first step towards the treat- 
ment of the complex dynamics of nucleation, thermal and 
photo fragmentations of PAH clusters in the interstellar 
medium. Our work could be extended by limiting the approx- 
imations done in the calculations of aggregation rates (impact 
parameter and initial energy distribution) and of evaporation 
rates (mainly the angular momentum). A next step could be 
to study larger clusters, especially those beyond the 400 car- 
bon atom limit. It would be important to determine whether 
such large clusters could survive inside a reflection nebula like 
NGC 7023. Unfortunately, the computational efforts for sim- 
ulating the evaporation of large clusters by molecular dynam- 
ics, to calculate their vibrational densities of states with Monte 
Carlo simulations, or to apply the IR/evaporation model will all 
tend to become extremely costly, at least one order of magni- 
tude higher than in the present study. Further approximations 
will then become unavoidable for a proper extension of the 
present work. 

Two important factors, namely ionization and density het- 
erogeneities, might enhance the stability of PAH clusters in 
PDRs. These factors have been qualitatively discussed here. 
Only singly-charged cationic clusters are expected to be more 
stable than neutrals. Conversely, multicharged species and an- 
ionic clusters should be less stable due to electrostatic repul- 
sion and electronic photodetachment processes, respectively. 
For the formation problem, it would be very useful to build a 
complete model to describe singly charged cationic assemblies 
of PAH molecules. 

To be chemically relevant, such a model would need to in- 
clude two main effects, that are not important for neutral clus- 
ters. Firstly, the charge will lead to l/r"* polarization contribu- 
tions, or even many-body implicit interactions between induced 
dipoles at higher order Secondly, the model should account for 
charge delocalization, the missing electron being able to hop 
between neighbouring molecules. In homogeneous molecular 
systems, the charge is not restricted to be carried by a single 
PAH, but may extend over se veral units, typically between two 
and four jPiuzzi et alJl2002l) . This phenomenon is known as 
charge resonance. Obviously, polarization and charge delocal- 
ization effects are not unrelated, but their many -body and quan- 
tum characters make them significantly harder to describe than 
the pairwise repulsion-dispersion-electrostatic interactions in 
the presently investigated neutral clusters. 

Thus, for cationic PAH clusters a more complex atomic 
model is necessary to take into account electron tranfer. A 
mo del along those lines was developed by Piuzzi et al. (2001 
and iBouvier et alJ (l2002h in order to investigate the structure 
of small clusters of aromatic molecules. The extension of this 
approach to large clusters and/or large PAH molecules will 
be computationally demanding, considering the extensive sam- 
pling required in the study of unimolecular evaporation. In the 
context of nucleation, the coupling between intermolecular and 
intramolecular modes will further burden the application of this 
model. With respect to neutral clusters, the theoretical model- 
ing of cationic PAH clusters is therefore expected to be more 
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involved. ElForts in this direction will be invaluable to improve 
our understanding of the physical and chemical evolution of 
photodissociation regions. 
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